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Abstract 

A numerical method to solve the fractional diffusion equation, which could also be 



< 

. easily extended to many other fractional dynamics equations, is considered. These 

fractional equations have been proposed in order to describe anomalous transport 
<*' I characterized by non-Markovian kinetics and the breakdown of Pick's law. In this 

paper we combine the forward time centered space (FTCS) method, well known 
for the numerical integration of ordinary diffusion equations, with the Griinwald- 

m 

' Letnikov definition of the fractional derivative operator to obtain an explicit frac- 

O ■ tional FTCS scheme for solving the fractional diffusion equation. The resulting 

r'"! ■ method is amenable to a stability analysis d la von Neumann. We show that the 

rS 

^ . analytical stability bounds are in excellent agreement with numerical tests. Com- 

parison between exact analytical solutions and numerical predictions are made. 
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1 Introduction 



Fractional differential equations have been a highly specialized and isolated 
field of mathematics for many years [1]. However, in the last decade there 
has been increasing interest in the description of physical and chemical pro- 
cesses by means of equations involving fractional derivatives and integrals. 
This mathematical technique has a broad potential range of applications [2]: 
relaxation in polymer systems, dynamics of protein molecules and the diffu- 
sion of contaminants in complex geological formations [3,4,5] are some of the 
most recently suggested [6]. 

Fractional kinetic equations have proved particularly useful in the context of 
anomalous slow diffusion (subdiffusion) [7] . Anomalous diffusion is character- 
ized by an asymptotic behavior of the mean square displacement of the form 



where 7 is the anomalous diffusion exponent. The process is usually referred 
to as subdiffusive when < 7 < 1. Ordinary (or Brownian) diffusion corre- 
sponds to 7 = 1 with Ki = D (the diffusion coefficient). From a continuous 
(macroscopic) point of view, the diffusion process is described by the diffu- 
sion equation ut{x,t) = D Uxx{x,t), where u{x,t) represents the probability 
density of finding a particle at x at time t, and where m^,^... is the partial deriva- 
tive with respect to the variables 77, C. . . From a microscopic point of view, the 
continuous description is known to be connected with a Markov process in 
which the microscopic particles (random walkers) perform stochastic jumps 
of finite mean and finite variance. In these conditions the central limit theo- 
rem holds for the sum of these jumps and Einstein's law for the mean square 
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displacement ensues [Eq. (1) with 7 = 1]. 

On the other hand, if an underlying non-Markovian microscopic process is 
assumed in which random walkers perform jumps at times chosen from a 
distribution with an algebraic long-time tail t"'''"^, then the diffusion process 
is anomalous [7,8]. In these circumstances the central limit theorem breaks 
down and one must apply the generalized Levy-Gnedenko statistics [7,9] which 
form the basis of Eq. (1). It turns out that the probability density function 
u{x, t) that describes these anomalous diffusive particles follows the fractional 
diffusion equation [7,10,11,12]: 

cP' 

-u{x, t) = qDI-^—u{x, t) (2) 

where oDI~^ is the fractional derivative defined through the Riemann-Liouville 
operator (see Sec. 2). Fractional subdiffusion-advection equations, and frac- 
tional Fokker-Planck equations have also been proposed [13,14,15,16] and even 
subdiffusion-limited reactions have been discussed within this framework [17]. 
In the mathematical literature, these equations are usually referred to as 
parabolic integro-differential equations with weakly singular kernels [18]. 

These current applications of fractional differential equations and many others 
that may well be devised in the near future make it imperative to search for 
methods of solution. Some exact analj^ical solutions for a few cases, although 
important, have been obtained by means of the Mellin transform [11,12] and 
the method of images [19]. The powerful method of separation of variables 
can also be applied to fractional equations in the same way as for the usual 
diffusion equations (an example is given in Sec. 4). Another route to solving 
fractional equations is through the integration of the product of the solu- 
tion of the corresponding non-fractional equation (the Brownian counterpart 
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obtained by setting 7 — > 1) and a one-sided Levy stable density [7,20,21]. 
However, as also for the Brownian case, the availabihty of numerical methods 
for solving (2) would be most desirable, especially for those cases where no 
analytical solution is available. One possibility was discussed recently by R. 
Gorenflo et al. [22,23,24] who presented a scheme to build discrete models of 
random walks suitable for the Monte Carlo simulation of random variables 
with a probability density governed by fractional diffusion equations. Another 
more standard approach is to build difference schemes of the type used for 
solving Volterra type integro-differential equations [18]. In this line, some im- 
plicit {backward Euler and Crank-Nicholson) methods have been proposed 
[18,25,26,27,28,29,30]. In this paper we shall use the forward Euler difference 
formula for the time derivative du/dt in Eq. (2) to build an explicit method 
that we will call the fractional Forward Time Centered Space (FTCS) method. 
For Brownian (7 = 1) diffusion equations, this explicit procedure is the sim- 
plest numerical methods workhorse [31,32]. However, for fractional diffusion 
equations, this explicit method has been overlooked perhaps because of the 
difficulty in finding the conditions under which the procedure is stable. This 
problem is solved here by means of an analysis of Fourier-von Neumman type. 

The plan of the paper is as follows. In Sec. 2 we give a short introduction to 
some results and definitions in fractional calculus. The numerical procedure 
to solve the fractional diffusion equation (2) by means of the exphcit FTCS 
method is given in Sec. 3. In this section we also discuss the stability and the 
truncating errors of the FTCS scheme. In Sec. 4 we compare exact analytical 
solutions with the numerical ones and check the reliability of the analytical 
stability condition. Some concluding remarks are given in Sec. 5. 
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2 Basic concepts of fractional calculus 

The notion of fractional calculus was anticipated by Leibniz, one of the founders 
of standard calculus, in a letter written in 1695 [1,7]. But it was in the next 
two centuries that this subject fully developed into a field of mathematics with 
work of Laplace, Cayley, Riemann, Liouville, and many others. 

There are two alternative definitions for the fractional derivative o-C^t"''^ of 
a function f{t) which coincide under relatively weak conditions. On the one 
hand, there is the Riemann-Liouville operator definition 

oD]-y(t) = J-tI- t dT , ^^^l , (3) 
" * ■'^ ^ r(7) dt Jo {t - t)i-t ' ^ ^ 

with < 7 < 1. For 7 = 1 one recovers the identity operator and for 7 = 
the ordinary first-order derivative. On the other hand, for any function f{t) 
that can be expressed in the form of a power series, the fractional derivative 
of order 1 — 7 at any point inside the convergence region of the power series 
can be written in the Griinwald-Letnikov form 

oDl-y{t) = Ihn ^ X: ^t'^fit - l^h), (4) 

where [t/h] means the integer part oit/h. The Griinwald-Letnikov definition is 
simply a generalization of the ordinary discretization formulas for integer order 
derivatives [1]. The Riemann-Liouville and the Griinwald-Letnikov approaches 
coincide under relatively weak conditions: if f{t) is continuous and f'{t) is 
integrable in the interval [0, t] then for every order < 1 — 7 < 1 both the 
Riemann-Liouville and the Griinwald-Letnikov derivatives exist and coincide 
for any time inside the interval [0, t] [1] . This theorem of fractional calculus 
assures the consistency of both definitions for most physical applications where 
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the functions are expected to be sufficiently smooth. 



The Griinwald-Letnikov definition is important for our purposes because it 
allows us to estimate qD]~^ f{t) numerically in a simple and efficient way: 



1 . ^ 

.D'rm = ^ E ^t'^fii - kh) + 0{hF) , (5) 



The order of the resulting approximation, p, depends on the choice of ui^ "'\ 
The approximation is of first order (p = 1) when u>^^ is the k-ih coefficient in 
the power series expansion of (1 — z)" [1,33], i.e., 

oo 

{l-zr = ^utz' (6) 

k=0 

SO that oj^f^ = (-l)*'^^) or, equivalently: 

4"' = 1, u,r = (l-^!^)4-\ 4=1,2,... (T) 

The coefficients ou^~"'^ of the second-order approximation (p = 2) can be 
obtained similarly [1,33]: 

These coefficients can be easily calculated using Fast Fourier Transforms [1]. 
However, for the fractional FTCS method discussed in this paper, we will show 
in the next section that nothing is gained by using second-order approxima- 
tions for the fractional derivative. Besides, the stability bound is smaller if we 
take the coefficients derived from Eq. (8). Finally, it is important to note that 
the error estimates given in (5) are valid only if either ^ 1 [1] or u{x,t) 
is sufficiently smooth at the time origin t = [34]. 
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3 Fractional ForwEird Time Centered Space method. 



We will use the customary notation xj — jAx, tm — mAt and u{xj,tm) = 
u^^^ ~ Uj^^ where C/j"*^ stands for the numerical estimate of the exact value of 
u{x, t) at the point {xj, tm)- In the usual FCTS method, the diffusion equation 
is replaced by a difference recurrence system for the quantities u^^^: 

with T{x,t) being the truncation term [31]. In the same way, the fractional 
equation is replaced by 

(m+l)_ M M_^ M irn) 

' ' ^K^ oD]-^ ./.^ ^ + T{x, t) . (10) 



At ^ " ' (Ax)2 

The estimate of the truncation term will be given in Sec. 3.2. Inserting the 
Griinwald-Letnikov definition of the fractional derivative given in Eq. (5) into 
Eq. (10), neglecting the truncation term, and rearraging the terms, we finally 
get the explicit FTCS difference scheme 



jj(m+l) ^ ^(m) +S^J2 4 - '^Uf ^ ^ + U]"^r^\ , (11) 

fe=0 

where = KjAt/[h^^^{Ax)'^]. In this scheme, Uj"''^^\ for every position j, 
is given explicitly in terms of all the previous states Uj^\ n = 0,1,..., m. 
Because the estimates C/j"*^ of u{xj,tm) are made at the times mAt, m — 
1,2,.. ., and because the evaluation of QDl~^u{xj, t) by means of (5) requires 
knowing u{xj, t) at the times nh, n — 0,1,2, . . ., it is natural to choose h — At. 
In this case, 

S, = A%(^. (12) 
The solution u{x, t) is a causal function of time with u{x, t) — ii t < 
(m^"'' = if n < —1), and we assume that the system is prepared in an 
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initial state ' — Uj . The iteration process described by Eq. (11) is easily 
implementable as a computer algorithm, but the resulting program is far more 
memory hungry than the elementary Markov diffusive analogue because, in 
evaluating C/j"*^^\ one has to save all the previous estimates Uj'^^^\ 
and Uj^i^^ forn = 0, 1, . . . m. However, the use of the short-memory principle 
[1] could alleviate this burden. Anyway, before tackling Eq. (11) seriously 
we must first discuss two fundamental questions concerning any integration 
algorithm: its stability and the magnitude of the errors committed by the 
replacement of the continuous equation by the discrete algorithm. 



3. 1 Stability of the fractional FTCS method 



We will make a von Neumann type stability analysis of the fractional FTCS 
difference scheme (11). We start by assuming a solution (a subdiffusion mode 
or eigenfunction) with the form = Cm^^'^''^^ where g is a real spatial wave 
number. Inserting this expression into (11) one gets 

^Cm- 45sin2 (^) f:4~'\m-k ■ (13) 

It is interesting to note that this equation is the discretized version of 

^?M = _4C.in=(«f£) „A'-^(t), (14) 

[with C — S{^t)'^] whose solution can be expressed in terms of the Mittag- 
Leffler function E^^—Xf^) [2,7]. This result is not unexpected because the 
subdiffusion modes of (2) decay as Mittag-Leffler functions [7] [e.g., see (30)]. 

The stabihty of the solution is determined by the behaviour of Cm- Unfortu- 
nately, solving Eq. (13) is much more difficult than solvin the corresponding 
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equation for the diffusive case. However, let us write 



Cm+l — 



(15) 



and let us assume for the moment that ^ = ^[q) is independent of time. 
Then Eq. (13) implies a closed equation for the amplification factor ^ of the 
subdiffusion mode: 



If 1^1 > 1 for some q, the temporal factor of the solution grows to infinity 
according to Eq. (15) and the mode is unstable. Considering the extreme 
value ^ = —1, we obtain from Eq. (16) the following stability bound on S^: 



The bound expressed in Eq. (17) depends on the number of iterations m. 
Nevertheless, this dependence is wak: for m > 1, approaches S^^^ = 
in the form of oscillations with small decaying amplitudes (see Fig. 1). Figure 
2, in which we plot A^" = 5*^2 ~ -S"^,! versus 7 for the first- and second-order 
coefficients, serves to gauge the amplitude of these oscillations. In fact, AS^ is 
the maximum value of S^^j^+i ~ ^^,mi ^ ^ 1 when the first-order coefficients 
(7) are used. We see that AS^ is certainly small for all 7. 

The value of hm^^oo S^,^ — can be deduced from Eq. (17) taking into 
account that the coefficients uj\^ are generated by the functions given in 
Eqs. (6) and (8). When the first-order coefficients given by (6) are used, one 
gets: 




(16) 




(17) 




1 1 



(18) 



2(1-0^-1^^-1 22-7 ■ 
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0.40 



CO 




m 



Fig. 1. First values of S^^m versus m for 7 = 1/2 when the first-order coefficients 
(circles) and second-order coefficients (squares) are used. The lines mark the corre- 
sponding limit values given by Eqs. (18) and (19) 

Similarly, when the second-order coefficients given by (8) are used, one gets: 



c-x 
•^7 



1 



2 (I - 2e + le) 



1-7 



43/2-7 ■ 



(19) 



We will verify numerically in Sec. 4 that the explicit integration method as 
given by Eq. (11) is stable when 



(20) 



sm^ (^) 

and unstable otherwise. As the maximum value of the square of the sine func- 
tion is bounded by 1, we can give a more conservative but simpler bound: the 
fractional FTCS method will be stable when 



{Ax) 



2—7 



(21) 
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Fig. 2. The difference ASj = S^2 ~ ^^i versus 7 when the first-order coefficients 
for Lo^j} "'^ [c.f. Eq. (7)] (solid line) and second-order coefficients [c.f. Eq. (8)] (dotted 
line) are used. 

The physical interpretation of this restriction is the same as for the diffusive 
case, namely, Eq. (21) means that the maximum allowed time step At is, up 
to a numerical factor, the (sub) diffusion time across a distance of length Ax 



Notice that the value of = 1/4^/^"'"' given by Eq. (19) is smaller than 
1/2^"'^ for any 7 < 1 (if 7 = 1 we recover the bound = 1/2 oi the usual ex- 
phcit FTCS method for the ordinary diffusion equation [31,32]). Consequently, 
the fractional FTCS method that uses a second-order approximation in the 
fractional derivative is "less robust" than the fractional FTCS method that 
uses the first-order coefficients a^^^"^^- Taking into account that the two meth- 
ods have the same precision (see Sec. 3.2) we note that nothing is gained by 



[c.f. Eq. (1)]. 
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using the fractional derivative with higher precision. Therefore, in practical 
applications, we will only use here the first-order coefficients (7). 



3.2 Truncating error of the fractional FTCS method 



The truncating error T(x, t) of the fractional FTCS difference scheme is [see 
(10)]: 



T{x,t)^ 



(m+l) _ (m) 



But 



and 



At 



(mJ-l') (m) 
Vj — tl ' 



(Ax)2 



' = «t + i^utt^t + 0{^tf 



(22) 



(23) 



1 m 

7-)i-7 L.M _ 9,,M I ..Ml _ \^ .,,1-7 



k=0 



+0 {hP) 
(24) 

so that, taking into account that u{x,t) is the exact solution of Eq. (2), we 
finally get from Eqs. (22), (23) and (24) the following result 

T{x, t) = 0{h^) + ^UuAt - oDl-^u,xxx + • • • (25) 



0{W) + 0(At) + 0{Axf . 



(26) 



Therefore, (i) assuming that the initial boundary data for u are consistent 
(as assumed for the usual FTCS method [31]) and (ii) assuming that u is 
sufficiently smooth at the origin t — Q [see remark below Eq. (8)], we conclude 
that the method discussed in this paper is unconditionally consistent for any 
order p because T(x, i) ^ as /i. At, Ax 0. As remarked above, in practical 
calculations is convenient to use h — /^t so that, due to the term 0(Ai) in 
(26), no improvements are achieved by considering higher orders than p = 1 
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in the fractional derivative. In is interesting to note that for the diffusion 
equation (7 = 1) it is possible to cancel out the last two terms in Eq. (25) with 
the choice At = {Axy/{6Kj), trhereby obtaining a scheme that is "second- 
order accurate" [31] . This is not possible for the fractional case because of the 
fractional operator. 



4 Numerical solutions and the stability bound on 5", 



The objective of this section is twofold: first we want to test the rehabihty 
of the numerical algorithm defined in Eq. (11) by applying it to two frac- 
tional problems with known exact solutions, and second we want to check the 
stability bounds obtained in Sec. 3.1. 



4-1 Numerical solution versus exact solution: two examples 



The fundamental solution of the subdiffusion equation in Eq. (2) corresponds 
to the problem defined in the unbounded space where the initial condition is 
u{x, t = 0) = d{x). This solution is called the propagator (or Green's function) 
and can be expressed in terms of Fox's H- function [7] : 



u{x, t) 



.jATrKyt-y 



0-10 
-^11 



\x\ 



(1-7/2,7/2) 



(0,1) 



(27) 



In our numerical solution we used the boundary conditions u{—L, t) — u{L, t) — 
with a sufficiently large L in order to avoid finite size effects. In Fig. 3 we 
compare the numerical integration results with the exact solution (27) for 
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X 

Fig. 3. Comparison between the exact subdiffusion propagator (lines) and the nu- 
merical integration results for 7 = 1/4 (squares), 7 = 1/2 (circles), 7 = 3/4 (trian- 
gles) and 7 = 1 (crosses) and t = 10. 

7 = 1/4, 1/2, 3/4, 1 at i = 10. The timestep used was At = 0.01 and 
Ax = ^K^{At)^/S^ with = I and = 0.28, 0.33, 0.4 and 0.5. All these 
values of are just below the stabihty bound (see Eq. (18)). The agree- 
ment is excellent except for 7 = 1/4 and x = 0, but this minor discrepancy is 
surely due to the large spatial cell Ax ~ 1.06 used in this case. 

We have also considered a problem with absorbing boundaries, m(0, t) — 
w(l, t) — 0, and initial condition u{x, t — — x{l — x). The exact analytical 
solution of Eq. (2) is easily found by the method of separation of variables: 
u{x,t) — X{x)T{t). We thus find Xn{x) = sin(n7ra;) and 

^^-K.XloD'rT, (28) 



14 



0.020 - 



0.015 



0.010 



0.005 



0.000 




Fig. 4. Numerical solution of the subdiffusion equation for the problem with absorb- 
ing boundary conditions, ^(0, t) = u{l, t) = 0, and initial condition u{x, 0) = x{l—x) 
versus the exact analytical result (lines) for t = 0.5. The solution u{x, t) is shown 
for 7 = 0.5 (triangles), 7 = 0.75 (squares) and 7 = 1 (circles). 

where A„ = nn, n— 1,2, The solution of Eq. (28) is found in terms of the 

Mittag-Leffler function [7]: 



(29) 



Imposing the initial condition we obtain 
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7r3 {2n + 1)3 



sin[(2n + l)7rx]E^[-K{2n + lynH'^] 



(30) 



In Fig. 4 we compare this exact solution with the results of the numerical 
integration scheme for 7 = 0.5, 7 = 0.75, and 7 = 1 for t = 0.5 and = 1. 
The values of used were = 0.33, 0.4, and 0.5 with Ax = 1/10, 1/20, 
and 1/50, respectively. The values of At for fixed and Ax stem from the 
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Fig. 5. Values of corresponding to the onset of instability versus the subdiffusion 
exponent 7. The solid line is the prediction of the Fourier-von Neumann analysis 
and the symbols denote the results of the numerical tests with the criterion in 
Eq. (32): stars, triangles and squares for the absorbing boundary problem with 
u{x, 0) = x{l — x) with M = 50, 100 and 1000, respectively, and circles for the 
propagator with M = 1000. 

definition of S^: 



Excellent agreement is observed for the three values of 7, it being slightly 
poorer for the smallest value which is not surprising because in this case the 
mesh size Ax = 1/10 used is the largest. 




(31) 
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4-2 Numerical check of the stability analysis 



We checked the stabihty bound on the value of the given in Eq. (18) in 
the following way. For a set of values of 7 in the interval [0, 1], and for values 
of starting at = 0.98S^ (in particular, for = 0.98/22-t + 0.001 n, 
n = 0, 1, 2, . . .) we apphed the fractional FTCS integration until step M. We 
say that the resulting integration for a given values of 7 and is unstable 
when the following condition is satisfied at any position j: 

^ ">S for any m = M - AM, M - AM + 1, . . . , M , (32) 



where S = 5 and AM = 10. This means that the numerical solution is consid- 
ered unstable if the quotient u^'^ /uj- becomes negative or larger than 2S at 
any of the last AM steps. (Of course, this criterion is arbitrary; however, the 
results do not change substantially for any other reasonable choice of 5 and 
AM.) Let S^'"" be the smallest value of = 0.98/2^-^ + 0.001 n that verifies 
the criterion (32). For the absorbing boundary problem we calculate these 
values using Ax = 1/2N with AT = 5 and M = 50, M = 100 and M = 1000. 
For the propagator, we calculate S^™ using M = 1000 and At = 5 x 10~^ in 
a lattice with absorbing frontiers placed at a; = —N^x and x = NAx with 
= 50. It is well known that for a lattice with 2N + 1 points (including the 
absorbing boundaries) the maximum value of sin (g Ax/2) in Eq. (17) occurs for 
qAx = {2N-1)ti/{2N), so that in Fig. (5) we plot Sf'^sm'^[{2N -1)ti/{4:N)]. 
We observe that for large M the stability bound predicted by Eq. (18) agrees 
with the result of the numerical test. The larger values obtained for smaller 
M mean that the method must be "very unstable" to fulfill our instability 
criterion in so few steps. The success of the numerical test is truly remark- 
able and supports the unorthodox application of the Fourier-von Neumann 
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-1.5 



-4 





X 



Fig. 6. The propagator u{x, t) for 7 = 1/2, = 1,5 = 0.36 and t = 0.005 (squares) 
and t = 0.05 (circles). The time step is At = 0.0005 and the spatial mesh Ax is 
obtained according to Eq. (31). The lines are plotted as a visual guide. 

stability analysis to the fractional FTCS scheme made in Sec. 3.1. 

In Fig. (6) we plot the numerical solution when = 0.36 > in the case 
of the propagator with 7 = 1/2. This kind of oscillatory behaviour in the 
unstable domain is typical for ordinary partial differential equations too. 



5 Concluding remcirks 

The availability of efficient numerical algorithms for the integration of frac- 
tional equations is important as these equations are becoming essential tools 
for the description of a wide range of systems [6] . In this paper we have dis- 
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cussed a numerical algorithm for the solution of the fractional (sub) diffusion 
equation (2). Although we have dealt with this particular equation, our pro- 
cedure could be extended to any fractional integro-differential equation by 
means of an obvious combination of the Griinwald-Letnikov definition of the 
fractional derivative [1,2,7] with standard discretization algorithms used in 
the context of ordinary partial differential equations [31]. Furthermore, the 
method (given its explicit nature) can be trivially extended to d-dimensional 
problems, which is not such an easy task when implicit methods are consid- 
ered. 

In our numerical method the state of the system at a given time t — mAt is 
given explicitly in terms of the previous states ai t — {m — I) At, . . . , At, by 
means of the FTCS scheme (11). We verified that for some standard initial 
conditions with exact analytical solution, namely, (a) the propagator in an 
unlimited system with u{x,t = 0) = 6{x) and (b) a system with absorbing 
boundaries and u{x,t = 0) = x{l — x), the present algorithm leads to nu- 
merical solutions which are in excellent agreement with the exact solutions. 
Using a Fourier-von Neumann technique we have provided the conditions for 
which the fractional FTCS method is stable. For example, if a first-order ap- 
proximation for the fractional derivative is considered, we have shown that 
the FTCS algorithm is stable if = K^{Aty /{Axf < 1/2^-^. For 7 = 1 the 
well-known bound S — DAt/{AxY < 1/2 ol the ordinary exphcit method for 
the diffusion equation is recovered. 

Concerning the implementation of the method we must remark that the evalu- 
ation of the state of the system at a given time step mAt requires information 
about all previous states aXt— {m — l)At, {m — 2)At, . . . , At, and not merely 
the immediately preceding one as occurs in ordinary diffusion. This is a conse- 
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quence of the non-Markovian nature of subdiffusion and implies the need for 
massive computer memory in order to store the evolution of the system, which 
is especially cumbersome in computations of long-time asymptotic behaviours. 
This could be paUiated by using the "short-memory" principle [1]. Another 
feature of the explicit numerical scheme is the interdependence of the temporal 
and spatial discrete steps for a fixed S^. If, as usual, one intends to integrate 
an equation with a given mesh Ax, then the corresponding step size At for a 
given < is of the order (AxY^"'. As a consequence. At could become ex- 
tremely small even for no too small values of Ax, especially when the problem 
is far from the diffusion limit, i.e., for small values of 7, so that the number of 
steps needed to reach even moderate times would become prohibitively large. 
In this case, the resort to implicit methods [18,25,26,27,28,29,30], stable for 
any value of At and Ax, is compulsory. 

This work has been supported by the Ministerio de Ciencia y Tecnologia 
(Spain) through Grant No. BFM2001-0718. 

References 

[1] I. Podlubny, Fractional Differential Equations, (Academic Press, San Diego, 
1999). 

[2] R. Hilfer, Ed., Applications of Fractional Calculus in Physics, (World Scientific, 
Singapore, 2000). 

[3] J. W. Kirchner, X. Feng and C. Neal, Nature 403 (2000) 524. 

[4] H. Scher, G. Margolin, R. Metzler and J. Klafter, Geophysical Research Letters 
29 (5) (2002) doi:10.1029/2001GL014123. 



20 



[5] B. Berkowitz, J. Klafter, R. Metzler and H. Scher, preprint 'cond-mat / 0202327| 
v2 (2002). 

[6] I. M. Sokolov, J. Klafter and A. Blumen, Physics Today 55 (11) (2002) 48. 
[7] R. Metzler, J. Klafter, Phys. Rep. 339 (2000) 1. 

[8] G. Rangarajan and M. Ding, Phys. Lett. A 273 (2000) 322; Phys. Rev. E 62 
(2000) 120. 

[9] B. H. Hughes, Random Walks and Random Environments, Volume 1: 
Random Walks (Oxford, Clarendon Press, 1995); Random Walks and Random 
Environments, Volume 2: Random Environments (Oxford, Clarendon Press, 
1995). 

[10] V. Balakrishnan, Physica A 132 (1985) 569. 

[11] W. Wyss, J. Math. Phys. 27 (1986) 2782. 

[12] W. R. Schneider and W. Wyss, J. Math. Phys. 30 (1989) 134. 

[13] R. Metzler, E. Barkai and J. Klafter, Phys. Rev. Lett. 82 (1999) 3563. 

[14] A. Compte, Phys. Rev. E 55 (1997) 6821. 

[15] A. Compte and M. O. Caceres, Phys. Rev. Lett. 81 (1998) 3140. 

[16] R. Metzler, J. Klafter and I. M. Sokolov, Phys. Rev. E 58 (1998) 1621. 

[17] S. B. Yuste and K. Lindenberg, Phys. Rev. Lett. 87 (2001) 118301; Chem. Phys. 
284 (2002) 169. 

[18] C. Chuanmiao and S. Tsimin, Finite Element Methods for Integrodifferential 
Equations (World Scientific, Singapore, 1998). 

[19] R. Metzler and J. Klafter, Physica A 278 (2000) 107. 

[20] E. Baraki and R. J. Silbey, J. Phys. Chem. B 104 (2000) 3866. 

21 



[21] I. M. Sokolov, preprint |cond-mat/0101232l (2001). 



[22] R. Gorenflo and F. Mainardi, Fract. Calculus Appl. Anal. 1 (1998) 167. 

[23] R. Gorenflo, G. De Frabritiis and F. Mainardi, Physica A 269 (1999) 79. 

[24] R. Gorenflo, F. Mainardi, D. Moretti, G. Pagnini and P. Paradisi, Chem. Phys. 
284 (2002) 521. 

[25] J. M. Sanz-Serna, SIAM J. Numer. Anal. 25 (2002) 319. 

[26] J. C. Lpez-Marcos, SIAM J. Numer. Anal. 27 (2002) 20. 

[27] CH. Lubich, I. H. Sloan, and V. Thome, Math. Comp. 65 (1996) 1. 

[28] W. McLean and V. Thome, J. Austral. Math. Soc. Ser. B 35 (1993) 23. 

[29] W. McLean, V. Thome and L. B. Wahlbin, J. Comp. Appl. Math. 69 (1996) 
49. 

[30] K. Adolfsson, M. Enelund and S. Larsson (preprint available at 
|http://www. math. chalmers.se/~stig/papers/adap.pdf I 



[31] K. W. Morton and D. F. Mayers, Numerical Solution of Partial Differential 
Equations, (Cambridge University Press, Cambridge, 1994). 

[32] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical 
Recipes in Fortran 77: The Art of Scientific Computing (second edition) 
(Cambridge University Press, Cambridge, 1992). 

[33] Ch. Lubich, SIAM J. Math. Anal. 17 (1986) 704. 

[34] R. Gorenflo, in: Fractals and Fractional Calculus in Continuum Mechanics, eds. 
A. Carpinteri and F. Mainardi (Springer Verlag, New York, 1997) p. 277. 



22 



